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We  present  a  new  numerical  method  for  calculating  an  evolving  2-D  Hele-Shaw  in¬ 
terface  when  surface  tension  effects  are  neglected.  In  the  case  where  the  flow  is  directed 
from  the  less  viscous  fluid  into  the  more  viscous  fluid,  the  motion  of  the  interface  is 
ill-posed;  small  deviations  in  the  initial  condition  will  produce  significant  changes  in 
the  ensuing  motion.  This  situation  is  disastrous  for  numerical  computation,  as  small 
round-off  errors  can  quickly  lead  to  large  inaccuracies  in  the  computed  solution.  Our 
method  of  computation  is  most  easily  formulated  using  a  conformal  map  from  the  fluid 
domain  into  a  unit  disk.  The  method  relies  on  analytically  continuing  the  initial  data 
and  equations  of  motion  into  the  region  exterior  to  the  disk,  where  the  evolution  prob¬ 
lem  becomes  well-posed.  The  equations  are  then  numerically  solved  in  the  extended 
domain.  The  presence  of  singularities  in  the  conformal  map  outside  of  the  disk  intro¬ 
duces  specific  structures  along  the  fluid  interface.  Our  method  can  explicitly  track  the 
location  of  isolated  pole  and  branch  point  singularities,  allowing  us  to  draw  connec¬ 
tions  between  the  development  of  interfacial  patterns  and  the  motion  of  singularities 
as  they  approach  the  unit  disk.  In  particular,  we  are  able  to  relate  physical  features 
such  as  finger  shape,  side-branch  formation,  and  competition  between  fingers  to  the 
nature  and  location  of  the  singularities.  The  usefulness  of  this  method  in  studying  the 
formation  of  topological  singularities  (self-intersections  of  the  interface)  is  also  pointed 
out. 
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Contract  No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications 
in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 
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1  Introduction 


The  displacement  of  a  viscous  fluid  by  a  less  viscous  fluid  in  a  Hele-Shaw  cell  has  been 
the  subject  of  intense  investigation  over  the  last  decade,  mainly  due  to  newly  discov¬ 
ered  mathematical  analogies  with  dendritic  crystal  growth,  directional  solidification, 
and  electro-chemical  growth.  The  original  motivation  behind  the  pioneering  work  of 
Saffman  &  Taylor  [36]  was  the  analogy  to  displacement  in  porous  medium.  Reviews 
by  Saffman  [38],  Bensimon,  Kadanoff,  Liang,  Shraiman,  &  Tang  [.5],  and  Homsy  [16] 
summarize  the  state  of  affairs  as  of  the  mid-eighties,  while  some  of  the  the  more  recent 
developments  are  reviewed  by  Pelce  [31],  Kessler,  Koplik,  &  Levine  [22],  Howison  [21], 
and  Tanveer  [42]  from  a  range  of  different  perspectives. 

In  this  paper,  we  shall  limit  our  investigations  to  channel  flow,  although  our  numer¬ 
ical  method  is  quite  general,  and  can  be  applied  to  other  geometries.  The  advantage 
of  this  restriction  is  that  our  description  can  be  specific.  Moreover,  the  phenomena 
displayed  in  channel  flow  have  been  studied  extensively,  and  are  representative  of  other 
Hele-Shaw  flows. 

A  steadily  advancing  flat  interface  in  a  channel  is  unstable  to  perturbations  when 
driven  by  the  less  viscous  fluid.  The  perturbations  grow  into  fingers,  but  the  subse¬ 
quent  behavior  depends  on  the  relative  strength  of  capillary  effects,  measured  by  the 
dimensionless  number  B  =  f{l2nVa^).  Here  a  is  the  surface  tension  coefficient,  b 
is  the  gap  width  of  the  cell,  a  is  the  channel  width,  V  is  the  speed  of  displacement  well 
in  front  of  the  interface,  and  /.i  is  the  viscosity  of  the  more  viscous  fluid  (the  viscosity  of 
the  less  viscous  fluid  is  assumed  negligible).  Numerical  computations  by  Tryggvasan 
&  Aref  [45,  46],  DeGregoria  and  Schwartz  [9],  Bensimon  [6],  and  Meiburg  and  Homsy 
[27]  show  competition  between  fingers  resulting  in  the  emergence  of  a  single  steady 
finger,  provided  B  is  greater  than  about  0.0004+  (but  not  greater  than  B^  =  0.025 
otherwise  the  interface  is  stabilized  by  capillary  effects).  For  still  smaller  B,  DeGre¬ 
goria  and  Schwartz  [9,  10]  and  Bensimon  [6]  find  that  the  finger  spontaneously  splits; 
this  can  be  induced  for  higher  values  of  B  by  introducing  perturbations  at  the  tip, 

^There  is  a  range  in  B,  depending  on  the  level  of  noise  in  an  experiment  or  numerical  calculation, 
in  which  a  transition  occurs  between  steadily  moving  fingers  and  continual  unsteady  motion.  A  typical 
range  is  0.0005  >  B  >  0.0002.  We  have  taken  0.0004  as  a  representative  value. 
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even  with  small  amplitude.  Experiments  also  reveal  the  emergence  of  tip  splitting  and 
side-branching  instabilities  (Park  &  Homsy  [28],  Tabeling,  Zocchi,  k  Libchaber  [41]). 
Bensimon  [6]  provides  numerical  evidence  supporting  heuristic  arguments  that  the  size 
of  the  perturbation  triggering  instability  decreases  with  B.  The  computations  of  Dai 
and  Shelley  [8]  (in  the  circular  geometry)  also  show  great  sensitivity  to  the  level  of 
numerical  precision  when  the  surface  tension  coefficient  is  small.  For  small  enough  5, 
even  noise  during  experiments  can  be  large  enough  to  set  off  a  pattern  of  continual 
tip-splitting  and  finger  competition  (  Maxworthy  [26],  Arneodo,  Couder,  Grasseau, 
Hakim  k  Rabaud  [2]).  Indeed,  when  capillary  effects  are  very  small,  it  appears  that 
the  pattern  is  fractal  (Maxworthy  [26],  Kopf-SiU  k  Homsy  [23],  Arneodo  et.  aL  [2]). 

Detailed  understanding  of  this  unsteady  behavior  is  limited.  In  particular,  the 
numerical  work  has  not  produced  a  clear  understanding  of  the  asymptotic  trends  as 
B  Q.  Unfortunately,  the  inclusion  of  surface  tension  in  the  usual  mathematical 
model  makes  theoretical  studies  very  difficult.  Instead  a  large  body  of  knowledge  has 
been  developed  for  the  initial-value  problem  when  5  —  0.  For  example,  Gustaffson 
[13,  14]  has  rigorously  proved  the  existence  of  a  solution  for  finite  time  starting  with 
analytic  initial  data.  Earlier,  Galin  [11]  and  Polubarinova-Kochina  [32]  considered 
the  mathematically  identical  problem  of  the  Darcy  model  for  ground  water  flow  and 
devised  analytical  techniques  to  obtain  exact  solutions  for  a  class  of  initial  conditions. 
These  were  apparently  well-known  in  the  Russian  literature  (see  Hohlov  [15]  and  How- 
ison  [21]).  Exact  solutions  due  to  Saffman  [37],  Howison  [18,  19,  20]  and  Shraiman 
k  Bensimon  [39]  can  be  seen  as  applications  of  these  techniques  though  these  results 
were  obtained  without  knowledge  of  the  earlier  Russian  work.  Howison  [21]  summa¬ 
rizes  the  relation  between  the  different  techniques.  Within  the  class  of  known  exact 
solutions,  there  are  finger  patterns  that  exist  for  all  times  and  exhibit  behavior  similar 
to  experimental  observation  (Patterson  [30]).  Further,  there  are  solutions  (Shraiman 
k  Bensimon  [39])  that  exist  only  for  a  finite  time  and  culminate  in  a  zero  angled  cusp 
at  the  interface.  Howison  [19]  uses  the  class  of  known  exact  solutions  to  point  out  that 
the  initial-value  problem  is  ill-posed;  it  is  possible  to  choose  an  initial  condition  for 
which  the  solution  exists  for  all  times,  whereas  there  is  a  neighboring  initial  condition 
for  which  the  interface  develops  a  cusp  after  a  finite  time. 
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In  essence,  these  theoretical  results  use  a  conformal  map  t)  which  maps  the 
interior  of  a  unit  semi-circle  in  the  C  plane  to  the  physical  flow  domain  of  a  channel 
(see  Figure  1).  The  location  of  the  free  surface  at  a  time  t  is  given  by  = 

x{0,  t)  -f  iy{6,  t)  for  ^  =  e‘^  on  the  arc  of  the  semi-circle.  The  equation  for  the  evolution 
of  (x{0,t),y{0,t))  results  from  the  usual  application  of  the  kinematic  and  dynamic 
conditions  at  the  interface.  The  symmetry  in  the  problem  allows  us  to  reflect  the 
solution  about  the  real  axis.  In  essence,  we  include  a  mirror  image  of  the  channel 
so  that  the  interface  may  be  considered  periodic.  Thus,  the  conformal  map  must  be 
analytic  inside  the  unit  disk,  aside  from  a  logarithmic  singularity  at  C  =  0,  but  it 
may  have  singularities  and  zeros  outside  it.  These  can  move  towards  the  boundary 
of  the  unit  disk,  |C|  =  1-  In  particular,  zeros  may  reach  the  boundary  in  finite  time, 
causing  a  cusp  to  form  on  the  interface.  The  origin  of  iU-posedness,  then,  is  that  small 
perturbations  can  introduce  a  zero  or  a  singularity  in  near  the  unit  disk,  which 
subsequently  moves  towards  it  and  reaches  it  quickly.  Following  the  work  of  Richardson 
[.35]  and  Lacey  [25],  Tan  veer  [43]  showed  that  all  singularities,  no  matter  what  type, 
will  move  towards  the  unit  disk,  while  preserving  their  type.  Indeed,  we  surmise 
that  round-off  error  in  traditional  numerical  calculations  introduces  singularities  at 
.some  distance  from  [C|  =  1;  these  subsequently  approach  the  unit  disk  and  lead  to 
the  random  pattern  of  tip-splitting,  side- branching  and  finger  competition  seen  in 
computations  for  increasingly  small  B. 

Ideally,  one  would  like  to  understand  the  Hele-Shaw  dynamics  for  a  small  non-zero 
B  perturbatively,  thus  exploiting  the  simpler  5  =  0  case.  There  are  several  hurdles 
in  accomplishing  this.  One  is  the  ill-posedness  of  the  5  =  0  problem  in  the  physical 
domain  |(,‘|  <  1.  Another  is  the  fact  that  information  in  the  5  =  0  problem  is  not 
complete.  We  do  not  have  an  exact  solution  to  the  general  initial  value  problem,  nor 
do  we  have  a  detailed  understanding  of  the  motion  of  singularities.  Furthermore,  it 
is  not  known  how  to  numerically  compute  solutions  to  the  5  =  0  problem  effectively. 
Conventional  numerical  simulation  in  the  physical  flow  domain  |(C|  <  1  suffers  from 
uncontrolled  growth  of  round-off  errors  (see  Aitchison  &  Howison  [1]).  Employing  a 
filtering  procedure  such  as  that  used  by  Krasny  [24]  in  the  Kelvin- Helmholtz  instability 
(another  ill-posed  problem)  allows  simulations  in  the  physical  domain  to  proceed.  Still, 
as  is  demonstrated  by  Dai  &  Shelley  [8],  the  choice  of  parameterization  has  a  strong 
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effect  on  the  accuracy  of  the  computation.  Moreover,  a  good  parameterization  is 
very  dependent  upon  the  initial  data.  For  data  with  a  general  collection  of  zeros  and 
singularities,  reliable  long  time  calculations  in  the  physical  domain  become  extremely 
difficult. 

In  a  recent  development.  Tan  veer  [43]  has  been  able  to  demonstrate  that  analyti¬ 
cally  extending  the  initial  value  problem  for  ^(4,  t)  into  the  region  exterior  to  the  unit 
disk  leads  to  a  well-posed  evolution  problem.  It  is  important  to  note  that  when  the 
initial  interface  location  is  known  only  to  a  finite  precision,  initial  data  0)  in  the 
extended  region  |4|  >  1  cannot  be  obtained  in  a  well-posed  manner.  In  effect,  the 
ill-posedness  of  the  dynamics  is  transfered  to  the  ill-posedness  of  extending  the  data 
from  Id  =  1  to  Id  >  1.  Nevertheless,  when  data  is  specified'm  |d  >  1  (say,  with  2(C»0) 
given  in  closed  form)  the  interface  evolves  without  sensitivity  to  initial  conditions. 

Tanveer’s  observation  addresses  the  first  hurdle  mentioned  above.  In  this  paper, 
we  address  the  second  hurdle  by  presenting  a  numerical  method  which  efficiently  solves 
the  initial  value  problem  in  the  expanded  domain  when  5  =  0.  Our  work  extends  the 
method,  developed  by  Baker  k  Tan  veer  [4],  to  include  the  trajectories  of  singularities 
of  the  form 

^  A{t)  iC  -  Csit))\  07^0,1,2,...  (1) 

explicitly  in  the  complex  plane  C  when  5  =  0.  Thus  we  are  able  to  assess  directly 
the  impact  of  the  close  approach  of  pole  and  branch  point  singularities  to  |d  =  1- 
We  do  find  that  singularities  induce  tip-splitting  and  side-branching,  and  that  the 
relative  strength  of  the  singularities  control  finger  competition.  Moreover,  our  studies 
draw  connections  between  the  parameters  specifying  the  singularities  and  the  result¬ 
ing  physical  behavior.  In  other  words,  for  given  initial  data  we  are  able  to  predict 
the  outcome  from  knowledge  about  the  initial  singularities  in  |(^|  >  1.  Within  this 
framework,  comparisons  with  experimently  observed  features  are  possible  by  studying 
a  random  ensemble  of  initial  conditions  in  |4“|  >  1,  subject  to  the  constraint  that  they 
describe  the  same  initial  interface,  up  to  some  ‘experimental’  error. 

Beside  the  explicit  treatment  of  singularities,  there  is  another  crucial  ingredient  in 
our  method.  An  analytic  function,  such  as  a  conformal  map,  is  determined  completely 
by  its  values  on  a  closed  curve.  Instead  of  advancing  the  conformal  map  on  |C|  =  1  in 
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time  as  in  standard  boundary  integral  methods,  our  method  advances  the  conformal 
map  on  a  much  larger  circle.  The  interpolation  of  the  conformal  map  to  |C|  =  1  is 
then  a  well-posed  operation.  In  essence,  knowledge  of  the  conformal  map  on  a  much 
larger  circle  corresponds  to  knowledge  of  the  conformal  map  on  |C|  =  1  to  a  higher 
degree  of  precision.  In  contrast,  extrapolation  of  the  conformal  map  to  a  larger  circle 
is  an  ill-posed  operation,  an  operation  not  to  be  attempted  numerically.  Since  our 
method  uses  the  discrete  Fourier  transform  to  evaluate  Laurent  series,  our  results  are 
spectrally  accurate  and  can  be  performed  in  0{N  log  N)  operations.  Furthermore, 
explicit  treatment  of  the  nearest  singularities  helps  avoid  point  crowding  typical  of 
conformal  maps. 

There  are  several  limitations  to  our  method.  First,  there  is  the  matter  of  the  initial 
data.  We  require  that  all  singularities  of  z,^  in  the  extended  domain  be  isolated  branch 
points  or  poles.  This  rules  out  some  types  of  initial  data;  for  example,  structures 
such  as  natural  boundaries  or  es.sential  singularities  are  not  allowed  in  the  extended 
domain.  Nevertheless,  the  kinds  of  data  we  do  allow  produce  a  wide  range  of  interfacial 
features  such  as  tip-splitting,  side-branching  and  finger  competition  which  are  similar 
to  experimental  observations. 

The  second  limitation  is  that,  at  present,  our  method  is  restricted  to  the  zero 
surface  tension  problem.  Given  this,  it  is  important  to  discuss  the  relationship  between 
B  =  0  solutions  and  those  for  B  >  0.  The  weU-posed  formulation  of  the  B  =  0 
problem  allows  the  effects  of  small  surface  tension  to  be  addressed  perturbatively.  It 
turns  out  that,  at  the  initial  time,  surface  tension  is  a  regular  perturbation  except  in 
the  neighborhood  of  zeros  and  certain  singularities  of  z^.  Tan  veer  [43]  has  examined 
the  effect  of  small  surface  tension  on  isolated  singularities  in  the  initial  data  of  the 
form  (1).^  Capillary  effects  on  singularities  with  a  <  —4/3  cause  only  a  regular 
perturbation.  Their  powers  are  unchanged  but  their  speed  and  strengths  A(t)  are 
modified  slightly  by  terms  proportional  to  B.  Singularities  with  —4/3  <  a  <  -1/2 
are  immediately  transformed  into  a  localized  cluster  of  —4/3  branch  point  singularities, 
though  the  behavior  (1)  is  still  relevant  in  an  outer  asymptotic  sense.  Thus,  for  times 
much  less  than  1/B  the  interface  behaves  as  though  it  is  unaffected  by  surface  tension, 

*Our  use  of  a  for  the  power  corresponds  to  Tanveer’s  [43]  choice  of  —0. 
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provided  these  singularities  do  not  come  within  a  distance  0{B^)  from  |4|  =  1  and 
zeros  in  are  far  from  the  unit  disk.§  For  initial  conditions  with  zeros  present,  recent 
evidence  (Siegel,  Tanveer,  &  Dai  [40])  shows  that  small  surface  tension  effects  can 
significantly  perturb  the  interface  in  0(1)  time.  This  can  happen  even  when  the  zeros 
do  not  impinge  on  the  unit  disk.  The  surface  tension  effects  occur  in  predictable  ways 
when  a  localized  cluster  of  —4/3  singularities  created  out  of  an  initial  zero  (termed 
daughter  singularity  by  Tanveer  [43])  come  within  an  neighborhood  of  the 

unit  disk.  However,  there  are  initial  conditions  for  which  daughter  singularity  effects 
do  not  occur  in  0(1)  time.  Our  computations  determine  the  leading  order  interfacial 
shapes  in  such  cases.  Perhaps  more  importantly,  the  method  presented  here  provides 
a  means  to  determine  if  and  when  the  zero  surface  tension  evolution  of  the  interface 
deviates  from  the  small  surface  tension  evolution. 

In  the  next  section,  we  describe  the  equations  upon  which  our  method  is  based. 
The  explicit  treatment  of  singularities  is  discussed  in  Section  3.  Then  we  describe  the 
numerical  method  in  Section  4.  In  Section  5,  we  present  tests  of  our  method,  and 
some  typical  results.  Our  conclusions  are  discussed  in  Section  6. 


2  The  Equations  of  Motion 

In  this  section  we  present  the  equations  which  govern  interfacial  flow  in  a  Hele-Shaw 
cell  without  the  details  of  their  derivation.  We  will  follow  closely  the  formulation  in 
Tanveer  [43].  Our  discussion  will  be  limited  to  flow  in  the  channel  geometry.  For  the 
equivalent  formulation  in  a  radial  geometry,  see  Tanveer  [43]. 

Consider  a  Hele-Shaw  cell  of  infinite  length  and  finite  width  a  in  which  air  of 
negligible  viscosity  is  pushing  a  viscous,  incompressible  liquid.  Introduce  the  conformal 
map  2(C,0  which  takes  the  interior  of  a  unit  semi-circle  in  the  C  plane  into  the  viscous 
fluid  region  of  the  channel,  which  lies  in  the  2  plane.  The  circular  arc  |(|  =  1  is  mapped 
to  the  interface,  and  the  diameter  is  mapped  to  the  channel  walls.  A  schematic  of  the 
mapping  is  provided  in  Figure  1.  Note  that  we  set  the  width  a  =  2. 

®When  singularities  are  initially  0(1)  distance  from  the  unit  disk,  they  can  move  to  within  0{B  2+«  ) 
in  0(ln  time. 
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The  functional  form  of  the  conformal  map  is  given  by 


= --In  c  +  *  + /(C,0  (2) 

TT 

where  /  is  analytic  in  an  open  set  which  contains  the  unit  semi-circle.  The  analyticity 
of  /  on  the  circular  arc  guarantees  the  smoothness  of  the  interface.  We  require  that 


Im{/}  =  0  (3) 

on  the  real  diameter  of  the  semi-circle  to  satisfy  the  condition  that  z  maps  this  diameter 
to  the  channel  walls.  In  addition,  we  assume  that  7^  0  in  a  region  containing  the 
unit  semi-circle.  The  Schwartz  reflection  principle  then  implies  that  /  is  analytic  and 
2^  7^  0  for  Id  <  1- 

The  fluid  velocity  u,  averaged  across  the  plate  gap,  satisfies  Darcy’s  law 


u  = 


where  /i  is  the  viscosity,  b  is  the  plate  gap,  and  p  is  the  pressure  (here  considered 
as  a  function  of  x  and  y).  Thus  {—b^ /12p)p  provides  a  velocity  potential  (j).  Incom¬ 
pressibility  implies  the  existence  of  a  stream  function  0.  We  can  therefore  introduce 
a  complex  potential  function  W{z,t)  =  <p{z,t)  -|-  ii>{z,t)  which  is  an  analytic  function 
of  2  in  the  fluid  region  of  the  channel.  Considered  as  a  function  of  C,  this  potential  is 
decomposed  as 

= --lnC  + f  +  w(C,<)  (4) 

TT 

where  uj  is  assumed  to  be  analytic  in  unit  semi-circle,  implying  (via  the  Schwartz 
reflection  principle)  its  analyticity  for  |d  <  1.  In  (4)  the  velocity  at  infinity  is  assumed 
to  be  1;  together  with  a  -  2,  this  choice  makes  our  variables  effectively  dimensionless. 
The  relation 


Im{a;}  =  0 


(5) 


is  required  to  hold  on  the  real  diameter  of  S,  and  is  the  mathematical  statement  of 
the  condition  that  there  is  no  flow  through  the  walls. 

The  interfacial  conditions  will  determine  the  evolution  of  the  map  z{(^,  t).  The  kine¬ 
matic  condition  states  that  the  normal  component  of  the  fluid  velocity  is  continuous 
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across  the  interface  and  implies  that  on  |C|  =  1 


Details  of  the  derivation  of  this  equation  are  available  in  Saffman  [37]  and  Richardson 
[35].  In  the  absence  of  capillary  effects,  the  dynamic  condition  requires  continuity  of 
pressure  across  the  interface.  When  the  viscosity  of  the  less  viscous  fluid  is  negligible, 


this  condition  gives 


Re{u;}  0 


on  1(1  =  1.  We  note  in  passing  that  a  more  complete  description  of  Hele-Shaw  flow 
introduces  more  complicated  interfacial  conditions  due  to  a  thin  film  left  behind  by  an 
advancing  interface.  However,  several  numerical  (Park  &  Homsy  [29],  Reinelt  [33,  34]) 
and  analytic  (Tanveer  [44])  studies  have  shown  that  many  salient  features  of  interfacial 
motion  described  by  the  more  detailed  model  are  captured  by  the  simpler  interfacial 
conditions  (6-7). 

From  (5)  and  (7)  it  is  immediately  apparent  that  u;  =  0  for  all  (  in  the  complex 
plane.  Thus,  (6)  becomes 


In  order  to  extend  (8)  into  the  region  |(|  >  1,  we  first  provide  an  analytic  extension 
of  the  conjugate  of  an  analytic  function  evaluated  on  (  (  1=  1.  Let 

oo 

nc)  =  E«nC” 

71=0 

be  an  analytic  function  in  |(|  <  1.  Then 

oo 

nc)  =  E«-c"  (9) 

71=0 

is  analytic  in  |C|  <  1,  where  the  overbar  denotes  complex  conjugation.  For  functions 
that  are  not  necessarily  analytic  at  the  origin,  though  analytic  on  a  segment  of  |(|  =  1, 
we  can  generalize  the  above  definition  of  F{()  through  the  relation  F{()  =  F{Q. 
F{l/()  is  then  the  analytic  extension  of  F  from  |  (  |—  1.  It  follows  that 


Re{F{0}  =  llF{0  +  Fil/0 


for  C  =  In  addition,  |F(C)P  =  F{QF(l/C)  on  C  =  e'®.  Thus  (8)  can  be  written  as 


\C^</  7r2^(C,t)ic(l/C,<) 

on  C  =  e'®- 

The  continuation  of  equation  (11)  into  the  domain  |C|  <  1  can  now  be  obtained  in 
a  straight  forward  manner  by  employing  the  Poisson  Integral  Formula.  In  particular, 
we  use  a  variant  of  the  standard  formula  which  gives  the  value  of  an  analytic  function 
in  the  domain  |C|  <  1  in  terms  of  its  real  part  evaluated  on  the  unit  circle.  Application 
of  the  formula  to  the  function  2</(C-^c)  yields  (with  appropriate  choice  of  imaginary 


constant) 


for  Id  <  1,  where 


f\c'\=i 


Equation  (12)  can  be  analytically  continued  into  the  domain  (d  >  1  by  deforming 
the  contour  in  the  usual  way,  producing  an  additional  term  from  the  residue  of  the 
pole  at  (.  Consequently,  we  have 

for  Id  >  1.  Note  that  /(C,  0  defines  different  analytic  functions  in  |d  <  1  and  |d  >  1. 

A  useful  alternative  form  for  (14)  may  be  obtained  by  employing  (10)  to  write  (8) 
as 

zt  ^  Cztjl/CJ)  2C  . 

Czc  ^c(i/C,0  ^c(i/C)' 

Since  5t(l/C)  and  i<(l/C)  are  analytic  in  |  C  |>  E  this  equation  provides  us  with  a 
second  expression  for 

CE(1/C,t) 


when  Id  >  1- 

In  order  to  make  the  structure  of  (12)  and  (14)  more  apparent,  define 

qi  = 

-2C 

ki^/cy 


Q2 
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(17) 

(18) 


Then  (12)  becomes 

o 

1! 

1 

(19) 

for  ]C]  <  1)  and  (14)  becomes 

1 

11 

(20) 

for  |(,‘(  >  1.  These  equations,  corresponding  to  equations  (3.1)  and  (3.3)  of  Tanveer 
[43],  have  the  advantage  over  (8)  of  allowing  studies  of  the  presence  and  influence  of 
singularities  of  in  If,-]  >  1.  Furthermore,  they  lead  to  the  development  of  a  new 
numerical  procedure  for  calculating  interfacial  motion  in  a  well-posed  fashion  for  a 
restricted  class  of  initial  conditions. 

Since  /((, 0  and  5^(1/C)  are  analytic  in  |C|  >  1,  so  too  are  and  q2{Cit), 

except  at  infinity,  where  qi  has  a  simple  pole,  i.e.,  grows  as  (.  Since  (20)  has  a 
form  analogous  to  a  first  order  hyperbolic  system  in  the  complex  plane,  (although  in 
reality  it  is  a  nonlinear  integro-differential  equation  with  coefficients  qi  and  q^  that 
depend  nonlocally  on  z),  the  analyticity  of  qi  and  q2  has  important  consequences  on 
the  presence  and  motion  of  singularities  of  z.  For  the  convenience  of  the  reader,  we 
provide  a  summary  of  what  is  known  (Richardson  [35],  Lacey  [25],  Tanveer  [43]): 

1.  There  is  no  spontaneous  generation  of  singularities  in  the  finite  complex  plane. 
Furthermore,  the  form  of  a  singularity  which  is  present  initially  in  the  region 
](^]  >  1  is  invariant  with  time.  Singularities  which  are  present  initially  at  infinity 
do  not  move  to  a  finite  location. 

2.  We  define  a  “characteristic”  in  the  plane  by 

^  =  (21) 

Let  (c{t)  =  Rcit)  then 

Tanveer  [43]  has  shown  that  the  right  hand  side  of  (22)  is  less  than  zero  when 
](^c|  >  1-  Consequently,  “information”  outside  of  the  unit  disc  flows  inward 
toward  )C|  =  1- 
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3.  This  is  particularly  true  for  isolated  singularities  of  the  form  (1).  The  location  of 
these  singularities  satisfies  (21),  i.e.  they  move  with  speed  — so  they 
move  towards  the  unit  circle.  Incidentally,  the  singularity  in  2:  which  is  present 
initially  at  =  0  does  not  move  in  time,  since  the  characteristic  speed  at  =  0 
(given  by  g'i(0,t))  is  zero. 

4.  Singularities  with  a  >  —1/2  reach  the  unit  circle  in  finite  time.  Singularities 
with  a  <  — 1/2  come  indefinitely  close  to,  but  never  reach  the  boundary  |4|  =  1. 

Several  properties  of  zeros  in  are  also  relevant.  For  example,  there  is  no  spon¬ 
taneous  generation  of  zeros  of  Z(^  in  the  finite  complex  plane,  although  zeros  which 
are  present  initially  at  infinity  can  move  to  a  finite  (  location.  When  a  zero  impinges 
on  Id  =  1,  it  produces  a  zero-angled  cusp  in  the  shape  of  the  interface;  see  Howison 
[19,  20]  and  Shraiman  &  Bensimon  [39]  for  exact  solutions  where  this  happens.  There 
is  no  physically  sensible  way  of  continuing  the  solution  in  time.  Unlike  singularities, 
the  zeros  of  zc^  do  not  move  generally  with  the  characteristic  speed  —qi  and  it  is  hard 
to  predict  if  a  given  zero  will  hit  the  unit  circle  or  stall  at  a  finite  distance  from  |d  =  1. 

In  this  paper,  our  interest  will  be  focussed  on  initial  conditions  containing  only 
isolated  singularities  of  form  (1)  with  a  <  -1/2,  so  they  do  not  reach  |d  =  1  in  finite 
time.  We  shall  also  pick  cases  where  zeros  in  2:^  do  not  reach  |d  =  1  during  the  time 
of  our  computations.  For  these  cases,  the  inclusion  of  surface  tension  acts  in  a  regular 
perturbative  manner,  so  the  results  indicate  what  can  be  expected  in  the  limit  of  weak 
surface  tension. 

3  Explicit  Treatment  of  Singularities 

Baker  &  Tanveer  [4]  use  (20)  to  solve  directly  for  z{(^t)  on  a  circle  in  the  (  plane  of 
radius  Rit)^  assuming  that  this  circle  does  not  contain  a  zero  of  z^  or  a  singularity  of 
f(^  t)  =  2;  —  ?!  4-  (2/7r)  Ind  The  advantage  of  computing  the  evolution  of  /((,  t)  on  the 
boundary  |d  =  R(t)  is  that  during  the  process  no  singularities  in  will  be  introduced 
in  Id  <  R{t).  Of  course,  singularities  in  may  be  present  outside  |d  =  R(i)t  bnt 
as  long  as  R{t)  shrinks  fast  enough,  singularities  will  not  intrude  into  |d  <  Rit)’-  the 
presence  of  zeros  must  be  checked  separately.  Furthermore,  a  function  which  is  analytic 
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in  ICI  </?(<)  can  be  determined  throughout  this  region  from  its  values  on  |C|  = 
through  evaluation  of  its  Taylor  series  -  this  step  is  well-posed  as  demonstrated  in  the 
next  section.  Thus,  the  interface,  2(|C|  =  l,t),  can  be  recovered  from  the  solution  on 

Id  =  R{t). 

Unfortunately,  R{t)  must  shrink  faster  than  the  approach  of  any  singularity  to 
Id  =  1.  Consequently,  the  computation  may  terminate  (when  R  =  1)  before  any 
interesting  structures  have  developed  on  the  interface.  For  initial  conditions  containing 
isolated  singularities  of  form  (1),  this  obstacle  may  be  removed  by  explicit  treatment  of 
the  singularities.  In  this  section,  we  derive  equations  in  which  singularities  in  |d  >  1 
are  treated  explicitly;  they  are  in  a  sense  ‘subtracted’  from  the  data  and  evolved 
separately.  In  the  next  section,  we  show  how  the  resulting  equations  can  be  solved  in 
an  efficient  and  well-posed  fashion. 

For  initial  data  containing  singularities  of  form  (1)  in  |d  >  I?  fke  solution  may  be 
written  in  the  form, 


Oj  -f  1 


2(c,o  = 


+  1- 


O'?  +1 


Qfj-f-l 


0(0. 


0(0. 


+  (?(C)  0  c 

w 


(23) 


where  |di  >  0  are  real  constants  (excluding  0,1,2,...);  and  £'00(^)’0)  ^ 
Here  d(0  J  =  0  •  •  -  ^  <^i  mark  the  location  of  singularities  on  the  real  line,  and  d  (0 
for  i  =  Ji  -h  1, . . . ,  </  indicate  singularities  with  non-zero  imaginary  part.  In  order  to 
satisfy  (3),  the  amplitudes  for  j  =  1,...,  Ji  must  be  real.  The  amplitudes 

E^(C,t)  for  i  =  Ji  +  1,...,  J  are,  in  general,  complex.  Conjugate  singularities  0(0 
with  amplitudes  E^C,t)  are  included  to  satisfy  condition  (3).  The  branch  cut  is 
chosen  so  that  the  argument  of  (1  —  C/O)  between  x  and  —ir:  we  have  chosen 
this  form  to  simplify  the  numerical  evaluation  of  the  interfacial  location  when  |C|  =  1- 
In  the  region  |C|  >  1,  there  can  be  other  singularities  at  0  which  are  not  explicitly 
represented.  However,  it  is  assumed  that  the  0  fhe  nearest  singularities  so  that, 
given  f,  the  inequality  |0|  >  R{t)  >  Kji  holds  for  all  j  =  1,...,J.  The  functions 
E^{C,t)  and  are  analytic  in  the  annulus  1  <  |C|  <  R(t).  If  all  singularities  are 

explicitly  represented,  these  functions  are  analytic  in  the  entire  region  exterior  to  the 
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unit  disk.  Unfortunately,  it  is  possible  for  and  G  to  have  singularities  in  |C|  <  1 
which  cancel  in  (23)  so  that  2:  +  (2/7r)lnC  is  analytic  in  |(|  <  1. 

In  order  for  (23)  to  be  a  solution  to  (20),  and  G  must  satisfy 

$  =  -9.(0(l).<)  12'1) 


-  =(«,  +  !)£  - - ^ 


Gt  -  q\  G(^- - -q-i  +  92-  (26) 

irQ 

In  the  case  of  initial  data  with  logarithmic  singularities  (poles  in  Z(^)  the  solution 
has  a  modified  form, 

.(C.o  =  in  (1  -  [.^(C.Oln  (1  - 

+  £^(C,/)lnfl-^)  +G(C,<)-^lnC  (27) 


where 

^  =  -qimU)  (28) 

EG  -  q,  E\  =  0  (29) 

f-,  r'  r,  \  ^  p)  f  ~  91(0(0.0  9i(0.0l  ...A, 

G.-9.G(=9.  +  g£  I - - (30) 

Of  course,  it  is  possible  to  have  a  combination  of  the  two  forms,  (23)  and  (27). 

Our  numerical  method  is  based  on  (24-26)  or  (28-30):  in  the  next  section,  we 
describe  how  the  solution  may  be  updated  in  a  well-posed  fashion. 


4  Numerical  Method 

Our  numerical  treatment  is  a  combination  of  tracking  the  positions  of  the  singular¬ 
ities  by  solving  (24)  (or  (28))  and  of  advancing  EG(,t)  and  G{(,t)  in  time  by  the 
method  of  lines.  Despite  appearances,  it  is  not  necessary  to  update  jU-’((,<),  and 
G((,t)  throughout  the  computational  domain,  an  annulus  1  <  |Ct  <  R{t)-  Knowl¬ 
edge  of  these  functions  on  the  circles  |C|  =  1  and  1C|  =  R{t)  is  sufficient  to  determine 
them  everywhere  inside  the  annulus.  Nevertheless,  we  must  use  a  special  procedure 
to  update  EH(^,t)  and  G{(,t)  in  a  numerically  stable  way. 
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We  introduce  a  decomposition  for  a  function  g{()  which  is  analytic  in  the  annulus 
by  writing  it  as  the  sum  of  a  function  which  is  analytic  inside  a  circle  ^  iZ 

and  a  function  g-{C)  which  is  analytic  outside  the  circle  |C|  =  1.  The  decomposition  is 
made  unique  by  requiring  (0)  =  0.  By  considering  the  Laurent  expansion  for  giC)^ 
we  note  that 


0 


9-iO  = 

k=—oo 

(31) 

9+iC)  =  '^SlkC^ 

(32) 

k=i 


where  the  Taylor  series  for  5+ (C)  will  have  a  radius  of  convergence  >  R,  and  the 
Laurent  series  for  g—{C)  converge  outside  |(^|  =  R-  !•  We  call  the  function 
inner  analytic  and  the  function  g^  outer  analytic. 

An  important  feature  of  an  inner  analytic  function  is  that  its  values  inside  a  closed 
curve  can  be  determined  from  the  values  on  the  curve  in  a  well-posed  fashion.  To  see 
this,  consider  an  inner  analytic  function  evaluated  on  a  circle  of  radius  R.  Let 

(  =  then 

oo 

!/+(«)  = 

k  =  l 

In  practice,  the  coefficients  gkR^  obtained  by  the  discrete  Fourier  Transform.  Since 
the  sum  converges,  these  coefficients  must  decay  with  increasing  k.  There  will  be  a 
value,  kr  say,  where  the  coefficients  reach  the  round-off  levels  of  a  computer.  If  suffi¬ 
cient  resolution  is  used  with  the  discrete  Fourier  Transform,  then  the  coefficients  with 
>  fcy.  will  contain  round-off  levels  and  not  their  actual  values.  This  has  important 
consequences  when  we  wish  to  evaluate  g^  on  a  different  circle  (^  =  re  .  Then 

oo  /  r\^ 

k=l 

so  each  coefficient  gkR''  is  multiplied  by  prior  to  using  the  discrete  Fourier 

Transform  again  to  evaluate  the  sum.  Ur  >  R,  then  round-off  errors  in  g^R'^  for  k  >  kr 
will  be  significantly  amplified,  even  to  the  point  that  all  accuracy  in  p+(C  =  re*^)  will 
be  lost.  However,  if  r  <  i?,  round-off  errors  are  decreased.  For  an  outer  analytic 
function  the  situation  is  reversed:  evaluation  on  r  >  iZ  (here  we  will  take  R  =  1)  is 
numerically  stable,  while  on  r  <  fZ  it  is  not.  The  exception  to  our  statements  is  when 
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there  is  a  known,  finite  number  of  terms  in  (31)  or  (32)  with  amplitudes  well  above 
round-off  levels.  Then  the  relative  errors  in  the  evaluation  of  or  g^  on  any  circle 
remains  small  when  we  use  only  the  finite  number  of  terms. 


4.1  Basic  Algorithm 

We  decompose  both  £•'((,  t)  and  G(C,  t)  into  inner  and  outer  analytic  parts  and  use 
the  method  of  lines  to  update  on  |C|  =  R-,  and  E^,  G-  on  |C|  =  1.  We 

obtain  evolution  equations  for  these  quantities  by  applying  the  projection  operators, 
H+f  f+  and  H-f  =  /_.  Thus,  from  (25)  and  (26),  we  have 


U-  {Ri{(,t)} 
H.LiEl  +  {a,  +  l} 


Ai(C,0-9i(0(0,0 

V  c  -  0 


0(0  )  j 


(33) 

(34) 


G-t  =  ~  0  +  92((i  o|  •  (3^) 

Before  giving  the  results  for  the  other  components,  we  note  that  qi  and  q2  have  very 
simple  decompositions.  From  (16)  and  (17): 


9i(Ci  0  —  C0(0  +  9i- 


(36) 


where  =  (qi{t)  contains  only  one  term,  whereas  (18)  shows  that  q2  is  outer  analytic 
only, 

q2{C,t)  =  <12- ■  (37) 

The  form  of  the  decompositions  (36)  and  (37)  is  very  beneficial  in  our  design  of  a 
well-posed  algorithm.  For  example,  we  use  the  facts  that  the  product  of  E!_^  with  qi 
and  the  product  of  Ei.  with  [91(C)  -  9i(Ca:)]/(C-  Cfc)  result  in  only  outer  analytic  parts 
to  simplify  the  'H+  projection  of  (25): 


=  n+{R2(c,t)} 

—  ^4- +  (rtj  +  1) 


( -  qi{Cj{t),t) 

V  c  -  0 


qiiCjG) 

0(0 


(38) 


The  resulting  equation  for  G+  simplifies  to 


G 


+  t 


(39) 


when  we  use  the  facts  that  -2qil{wC),  92  and  91  are  outer  analytic. 
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The  following  properties  can  be  easily  deduced  from  (34-35)  and  (38-39): 

1.  The  evolution  of  does  not  depend  explicitly  on  EL-  As  a  result,  if  £’•^((’,0) 
satisfies  £^((,0)  =  0,  then  £^(C,0  “  0  for  all  t.  The  same  result  holds  for 

G+(C,0. 

2.  If  £:5_(C,0)  contains  a  finite  number  of  modes  in  (32)  with  highest  wavenumber 
k  —  M,  then  no  modes  with  wavenumber  k  >  M  will  be  generated  in  E^{(,t), 
The  same  result  holds  for  G4.. 

Similar  results  hold  for  the  decomposition  of  E^  and  G  in  (29)  and  (30)  for  loga¬ 
rithmic  singularities. 

To  apply  the  method  of  lines  to  (34-35)  and  (38-39),  we  assume  that  we  know  £^ 
and  at  N  evenly  spaced  points  on  the  circle  r  =  £(/),  and  Ei_  and  G^  at  N  evenly 
spaced  points  on  the  unit  circle.  Actually,  it  is  not  necessary  to  use  the  same  number 
of  points  on  both  circles,  but  we  make  that  assumption  for  ease  of  presentation.  First, 
we  describe  how  to  evaluate  the  right  hand  side  of  (34)  and  (35)  on  the  unit  circle.  For 
the  moment  we  assume  that  and  q2  are  known  there;  we  describe  their  computation 
in  detail  in  the  next  sub-section. 

We  compute  the  coefficients  CkR^  in  the  representation 

N/2 

k=i 

by  use  of  the  discrete  Fourier  Transform.  The  coefficients  are  obtained  by  a  division 
by  R^  and  then  used  to  evaluate  E^  and  E^^  on  the  unit  circle  by  use  of  the  discrete 
Fourier  Transform,  where 

N/2 

(40) 

k=l 

The  function  EL^  is  also  determined  on  |C|  =  1  by  using  the  discrete  fourier  Transform. 
Then  we  have  aU  the  information  needed  to  evaluate  £1  in  (34)  at  N  evenly  spaced 
points  on  the  unit  circle.  To  execute  the  projection  _ ,  we  use  the  discrete  Fourier 
Transform  to  calculate  the  coefficients  in  the  representation 

N/2 

k=-N/2 
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Then  we  set  to  zero  all  the  coefficients  with  k  >  1.  The  result  is  an  outer  analytic 
function,  which  is  the  forcing  term  in  (34).  Balance  of  the  Fourier  coefficients  of  like 
modes  on  the  left  and  right  hand  sides  of  (34)  then  yields  a  set  of  equations 


dEi 


dt 


-(0  =  h{t) 


for  the  Fourier  coefficients  El[t)  [k  <  0)  of  Ei,  In  an  equivalent  procedure,  we  may 
use  G_|_,  known  on  r  =  jR(/),  to  determine  and  use  G_,  known  on  the  unit  circle, 
to  get  Then  the  forcing  term  for  (35)  is  computed  by  executing  the  projection 
7i-  as  described  above. 

From  our  discussion  on  the  importance  of  “characteristics”  defined  by  (21),  we 
anticipate  that  E^^  and  G+  must  be  evaluated  on  a  circle  with  a  radius  that  is  collapsing 
at  rate, 


R  dt 


max  Re 
K\=R(t) 


(41) 


Then,  information  outside  the  circle  |(|  =  R{t)  will  not  cross  the  boundary  into  the 
annulus,  1  <  |C|  <  R(t).  We  include  an  advection  term  that  accounts  for  the  change 
in  E].  due  to  the  change  in  R{t):  thus  (38)  becomes 


dE\ 

dt 


=  {«.«,<)  + If  4<} 


(42) 


A  similar  term  is  needed  in  the  equation  (39)  for  G'+. 

We  will  describe  in  detail  in  the  next  sub-section  how  qi  and  92  may  be  evaluated 
on  Id  =  R{t}.  The  quantities  jFI.,  can  be  evaluated  by  spectral  techniques  as 
described  in  (40).  The  expression  contained  within  brackets  on  the  right  hand  side 
of  (42)  can  then  be  evaluated.  The  Fourier  coefficients  of  the  right  hand  side  of  (42) 
are  obtained  using  the  discrete  Fourier  Transform,  and  the  projection  H+  is  done  by 
zeroing  all  modes  with  k  <  0.  Equating  like  Fourier  modes  then  leads  to  a  set  of 
evolution  equations  for  El(t),  the  Fourier  coefficients  of  for  k  >  0.  The  same 
procedure  may  be  used  on  the  equation  for  G+. 

Any  suitable  ODE  Solver  may  be  applied  to  the  evolution  equations  for  Q,  E^,  and 
Gk-  For  the  results  used  in  this  paper,  we  use  the  standard  fourth-order  Runge-Kutta 
method  with  fixed  step-size. 

Since  z((,  t)  -f  (2/5r)ln  C  has  a  Taylor  series  expansion  about  C  =  0,  all  the  negative 
terms  in  the  Laurent  expansion  for  G(C,<)  must  cancel  all  the  negative  terms  in  the 
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Laurent  expansions  of  the  sums  in  (23).  This  will  be  only  approximately  true  in  the 
discrete  calculation,  so  it  can  be  used  as  a  check  for  numerical  accuracy.  Although  it 
is  not  necessary  to  compute  to  update  the  interface,  we  have  done  so  in  order 

to  employ  this  check. 

For  certain  initial  conditions,  our  algorithm  is  greatly  simplified.  For  example, 
when  E^.  has  only  a  few  terms  in  its  Taylor  series  expansion,  then  it  is  possible  to 
write  down  by  hand  the  evolution  equations  for  the  Fourier  coefficients,  since  the 
products 

C  ^ 

R  dt 

will  contribute  only  a  few  terms  due  to  the  simple  decomposition  (36)  for  qi-  In 
particular,  if  E\.,  Ei,  G+,  and  G-  are  time-dependent  constants  in  C,  then  only  Ei. 
need  be  computed. 

4.2  Computation  of  qi  and  q2 

Here  we  provide  details  of  the  computation  of  qi  and  q2  in  |C|  >  1  by  using  (13),  (17), 
and  (18).  First,  we  calculate  E].{C,t),  E[^{<^,t),  and  Ei^{CJ)  on  |C|  =  1  by  means  of 
the  discrete  Fourier  Transform  as  described  in  the  previous  subsection.  These  values 
are  used  to  compute  l/z,^(Q  =  pooiO/D{Q  on  |C|  =  1  where 

m  =  -  Ew(c)  [(»,  + 1)^  -  £|(c)  (i  - 

jWi-l-l  L  \  ^3/ \ 

+  GciO  PooiO  (43) 

and 

Prsio = n  (i  -  z)  n 

This  form  for  l/z^{C)  is  obtained  by  differentiating  (23)  with  respect  to  (,  then  fac¬ 
toring  out  l/poo(C)  taking  the  reciprocal.  We  have  found  it  necessary  to  compute 
l/z^(Q  using  this  expression  to  obtain  accurate  results  for  qi  when  singularities  are 
very  close  to  |C|  =  1.  The  function  1/5^(1/C)  can  be  computed  on  |C|  =  1  by  taking 
the  conjugate  of  1/^(;(C)-  Consequently,  the  function  q2  can  be  computed  on  the  unit 
circle  using  (18). 


9i(C,0 

c- 


9i(Cji  f ) 
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The  function  q\  is  computed  as  follows.  We  write 


c\  OO 

=  +  (44, 

on  (  =  We  determine  approximations  to  the  N  coefficients  |di, . . . ,  dyv/2}? 

|di, . . . ,  dyv/2}  by  means  of  the  discrete  Fourier  Transform.  Upon  substitution  of 
(44)  into  (13)  and  computation  of  the  residues,  one  finds  that 


Qi  ^  -C 


N/2 

do  +  2  ^  dk  C  ^ 

Ar=l 


(45) 


Note  that  dj  is  real  (i.e.,  dj  =  dj)  for  the  channel  geometry,  as  a  consequence  of  ^  being 
real  along  the  real  axis  of  The  functions  and  92  can  be  analytically  continued  out 
to  the  circle  \(\  =  R{t)  by  the  ideas  expressed  in  (40).  Similarly,  q\  can  be  evaluated 
at  the  location  (j{t)  of  the  singularities. 


4.3  Computations  for  Extremely  Close  Singularities 

Many  of  the  interesting  interfacial  features  revealed  by  our  calculations  occur  when 
singularities  approach  the  unit  circle  within  a  distance  much  less  than  machine  preci- 

■i 

sion.  To  track  these  singularities  reliably  when  they  are  that  close  to  the  unit  circle 
we  must  express  their  location  as 

0  =  (1 + 

The  time  evolution  of  the  quantities  Oj{t)  and  Aj(/)  ==  1/Sj  are  then  determined 
numerically  using  the  equations 

We  use  Aj  instead  of  Sj  so  that  singularities  can  be  allowed  to  come  extremely  close 
without  the  worry  of  them  jumping  inside  the  unit  circle  due  to  time-stepping  errors. 

The  computation  of  q\{(j^t)  is  more  delicate  when  Q  is  extremely  close  to  the  unit 
circle,  since  the  quantities 

1  _  1 
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in  (45)  will  be  approximated  as  1,  leading  to  inaccuracies  in  the  computation.  In  order 
to  obtain  a  better  value  for  1/Cj  we  write 


and  compute  defined  by 


1  +  fki^j) , 


using  the  recursion 


MSi)  -  1  +  .: . -■ 


The  substitution  of  (46)  into  (45)  then  yields  qi{Cj,t)  =  9l(e'^^0  +  where 


Then, 


+  2  X]  4  [^j(l  +  fk)  +  fk]  e' 

\  kz=zl 


where  we  have  used 


qiiCjy 


I  C  J  ^^dOki^/O 

when  Id  =  1,  which  follows  from  (13)  and  (17).  By  using  (48),  we  avoid  any  problems 
with  round-off  errors.  On  the  other  hand,  Im{gi/d}  is  computed  with  qi  directly 
determined  from  (47). 


5  Numerical  Results 

We  use  a  known  exact  solution  due  to  Saffman  [37]  to  help  validate  our  method: 
z{C,t)  =  i  +  d{t)-^  In  C+i  In 

where  1  <  «(0)  <  oo.  The  functions  d{t)  and  a{t)  are  determined  by 


ait)  =  1  +  Kie' 


d{t)  =  A'o  -b  2<  -I-  —  In  [1  -I-  /v'le" 
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where  Kq  and  A’l  are  constants  which  are  determined  by  the  initial  conditions.  From 
(49),  it  is  clear  that  has  poles  at  ^1,2  =  zeros.  The  form  of  the  initial 

data  (27)  used  in  our  numerical  method  is  made  to  correspond  to  (49)  by  setting 
£i((^,0)  =  £^(C,  0)  =  I/tt  and  G(C,0)  =  i  +  d(0).  We  pick  d(0)  so  that  the  initial 
profile  has  zero  mean  height.  The  functions  and  G{(^t)  are  constants  in  (  for 

this  problem,  and  therefore  only  a  single  Fourier  mode  (namely  the  constant  mode) 
is  necessary  to  specify  them.  Nevertheless,  we  use  512  modes  in  our  computation  in 
order  to  check  the  numerical  stability  of  our  algorithm  for  computing  and 

G(C,0« 

The  results  of  our  numerical  method  for  a(0)  =  2.0  are  shown  in  Figure  2.  We  use 
a  time  step  of  At  =  0.005.  At  ^  =  3,  the  calculation  gives  the  positions  of  the  poles  as 
(1^2  =  ±1.0000000244216,  whereas  the  exact  positions  are  (1^2  =  ±1.0000000244215. 
Although  the  difference  is  less  than  10“^^,  the  exact  and  the  numerically  computed 
profiles  differ  by  only  a  little  less  than  10^^.  However,  we  find  no  growth  of  round-off 
errors  in  the  modes  comprising  E\(^t)  and  In  this  example  we  set  R{0)  = 

2,000,  and  stop  the  calculation  at  i  =  3  when  R  is  nearly  1.  Because  there  are  no 
zeros  in  Z(^  we  can  take  i?(0)  much  larger  and  run  for  even  longer  times. 

There  are  no  known  exact  solutions  for  channel  geometry  when  branch  point  sin¬ 
gularities  are  located  in  the  region  \(\  >  1  (although  exact  solutions  with  branch  point 
singularities  have  been  found  for  sector  geometry;  see  Tu  [47]).  We  replace  the  two 
pole  singularities  in  Z(^  of  (49)  by  two  branch  point  singularities  of  equal  power  and 
amplitude,  and  locate  them  symmetrically  on  the  real  axis  of  Our  choice  is  mo¬ 
tivated  by  the  desire  to  understand  the  role  of  the  power  of  a  singularity  in  on 
the  shape  of  the  interface.  Thus  we  select  initial  data  corresponding  to  (23)  with 
Ji  =  2,  oi  =  02,  and  Ci(0)  =  “(2(0).  In  Figure  3,  oi  =  -4/5,  and  in  Figure  4, 
oj  =  —4/3.  Initial  values  for  the  singularity  positions  Cj{0)  and  amplitudes 
are  determined  through  experimentation;  those  that  lead  to  solutions  in  which  the 
interface  becomes  well-deformed  before  zeros  impinge  on  the  unit  disc  are  selected 
for  presentation.  In  this  manner,  the  initial  singularity  positions  and  amplitudes 
are  chosen  as  Ci(0)  =  — C2(0)  =  1*2,  E^{(^0)  =  E^{(^0)  =  1.94  for  Figure  3,  and 
Ci(0)  =  -(2(0)  =  3.3,  £i(C,0)  =  E^{C,0)  =  -2.2  for  Figure  4.  In  these  and  all 
subsequent  computations  we  set  G(C,0)  =  i  +  c^  where  c  is  a  real  constant  selected  to 
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produce  an  initial  profile  with  zero  mean  height.  Recall  that  for  constant  initial  data 
such  as  this,  only  negative  Fourier  modes  in  and  G  are  generated.  Consequently, 
we  only  need  to  solve  equations  (34)  and  (35).  The  positive  Fourier  modes  are  set  to 
zero  after  each  time  step  in  the  calculation.  We  use  512  points  on  the  unit  circle  in 
the  C  plane,  and  a  time  step  of  A<  =  0.0005. 

We  check  our  numerical  results  by  comparing  them  with  the  results  obtained  when 
none  of  the  singularities  outside  the  unit  disk  are  represented  explicitly.  To  perform 
the  latter  calculation,  we  use  the  code  to  solve  for  6'+(4,  i)  +  (jo(0  ^  circle  of  radius 

R{t)  that  is  closer  than  the  smallest  \Cj{t)\.  Consequently,  we  shrink  the  radius  R(t) 
as  described  in  section  4.  While  the  calculation  must  be  stopped  before  long  because 
R  <  1,  the  output  provides  a  useful  check  on  the  present  results  for  some  period  of 
time.  In  addition,  we  decrease  the  time  step  and  increase  the  number  of  modes  until 
there  are  no  detectable  differences  in  the  solution  within  plotting  accuracy. 

Unlike  the  previous  example,  we  cannot  be  sure  that  zeros  do  not  approach  the 
unit  disk.  So  we  monitor  the  presence  of  zeros  by  computing  the  integral 


1 


(50) 


which,  according  to  the  argument  principle,  equals  the  number  of  zeroes  minus  the 
number  of  poles  of  inside  the  unit  disk.  As  long  as  there  are  no  zeros  of  in 
|(^|  <  1,  this  integral  will  equal  -1,  due  to  the  simple  pole  of  at  C  =  0.  For 
the  results  presented  in  Figures  3  and  4,  the  value  of  remains  within  1%  of  -1 
throughout  the  length  of  the  computation.  Note  further  that  if  a  zero  in  z^  is  located 
near  the  unit  disk,  the  close  presence  of  a  pole  singularity  in  the  integrand  of  (50) 
causes  a  loss  of  accuracy  in  its  numerical  evaluation.  Thus  we  feel  confident  that  no 
zeros  are  very  near  to  the  unit  disk  up  until  the  times  of  the  final  profiles  in  Figures 
3,4,  and  5. 

The  contrast  between  Figures  3  and  4  is  quite  striking.  The  higher  value  of 
a  —  —4/b  in  Figure  3  correspond  to  a  singularity  of  weaker  effect  in  that  the  bulge  of 
fluid  at  the  base  of  the  finger  is  less  rounded,  causing  a  finger  that  is  more  pointed. 
The  lower  value  of  o  =  -4/3  in  Figure  4  produces  a  more  spherically  shaped  bulge  of 
fluid,  causing  a  thinner  neck  at  the  base  of  the  finger.  In  general,  our  experiments  with 
various  powers  a  show  that  for  powers  a  >  -1,  more  pointed  fingers  are  produced 
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from  the  less  pronounced  bulges  of  fluid  at  their  bases.  For  a  <  —  1,  the  fingers  have 
parallel  sides,  but  their  bases  have  thinner  necks  as  a  consequence  of  more  spherically 
rounded  bulges  of  fluid  there. 

A  representative  Fourier  spectrum  from  the  calculation  with  a  =  —4/3  is  pre¬ 
sented  in  Figure  5.  There  is  no  sign  of  spurious  growth  in  the  large  |A;|  modes  or  any 
other  indications  of  numerical  instability,  although  it  is  necessary  to  evaluate  q\  using 
expression  (43)  in  order  to  avoid  such  growth.  The  rise  in  the  tail  of  the  spectrum 
with  respect  to  time  is  due  to  singularities  in  the  functions  E\  G,  gi,  and  q2  at  po¬ 
sitions  C  =  1/Ci  and  C  =  1/(2  inside  the  unit  disc.  Since  ^  -  (2/7r)ln(  is  analytic  in 
Id  <  1,  the  singularities  in  and  G  must  cancel  out  in  the  expression  for  z{Q,  (23). 
Nevertheless,  these  singularities  do  affect  the  numerical  computation  of  the  quantities 
E^,  G,  qi,  or  q2.  Since  the  singularities  move  toward  (  =  ±1  from  inside  the  unit 
disc  as  Cl, 2  ^  ±1  from  the  outside,  a  large  number  of  Fourier  modes  are  required 
to  obtain  accurate  representations  for  these  quantities  when  the  singularities  are  very 
close  to  the  unit  circle.  Another  consequence  of  the  close  approach  of  the  singularities 
to  the  unit  circle  is  the  need  for  much  smaller  time  steps  to  maintain  accuracy.  This 
is  also  true  when  zeros  in  2:^  get  close  to  the  unit  disk.  In  general,  it  is  either  the 
close  approach  of  singularities  in  E^  and  G  or  zeros  in  to  |d  =  1  that  force  us  to 
terminate  our  calculations. 

If  we  ignore  the  loss  of  accuracy,  we  can  continue  the  calculations  shown  in  Figures 
3  and  4  a  little  further  in  time.  We  see  evidence  that  zeros  are  impinging  on  the  unit 
disk.  For  the  calculations  associated  with  Figure  3,  a  zero  is  approaching  |d  =  1  at 
a  point  corresponding  to  the  tip  of  the  finger.  Thus,  we  expect  the  finger  to  form  a 
cusp  in  finite  time.  For  Figure  4,  two  zeros  approach  |d  =  1  at  points  corresponding 
to  the  tops  of  either  side  of  the  finger,  so  we  expect  cusps  to  form  at  these  positions. 
Asymptotic  theory  suggests  that  the  initial  zeros  will  give  rise  to  localized  clusters 
of  -4/3  singularities  (daughter  singularities)  when  0  <  jB  <<  1.  The  leading  order 
motion  of  each  of  these  clusters  satisfies  equation  (21),  i.e.,  a  cluster  located  at  Q 
moves  with  speed  “Q'i(Q(0’0-  If  such  a  cluster  comes  close  to  |d  =  1,  it  can  cause 
the  interface  to  deviate  significantly  from  the  B  ~  0  solution. 

We  briefly  consider  the  influence  of  non-zero  surface  tension  on  interfacial  shapes 
by  comparing  the  zero  surface  tension  solution  to  that  for  small  surface  tension.  The 


23 


non-zero  B  solution  is  obtained  using  the  boundary  integral  method  developed  by  Hou, 
Lowengrub  k  Shelley  [17].  In  the  first  case  (Figure  6a),  we  consider  initial  data  with 
zeros  of  initially  placed  at  infinity,  and  with  singularities  which  satisfy  aj  <  —4/3.^ 
The  value  of  surface  tension  is  set  to  B  =  .0025.  From  asymptotic  theory,  it  is  expected 
that  the  addition  of  a  small  amount  of  surface  tension  will  make  little  difference  in  the 
evolution  of  the  interface  for  at  least  0(ln  B)  time.  The  actual  agreement  between  the 
interfacial  shapes  is  quite  remarkable:  the  two  solutions  are  indistinguishable  over  the 
entire  length  of  the  run.  The  agreement  is  unaffected  by  using  real  singularities  which 
only  satisfy  aj  <  -1/2,  since  these  singularities  behave  as  though  they  are  essentially 
unaffected  by  capillary  effects  for  the  time  of  the  computation. 

In  contrast,  when  the  initial  data  of  Figure  4  is  used,  the  difference  between  the 
B  =  0  and  B  =  .0025  shapes  is  very  significant.  As  seen  in  Figure  6b,  the  B  >  0 
finger  eventually  diverges  from  the  corresponding  zero  surface  tension  solution  and 
approaches  a  broad,  steadily  propagating  finger.  The  broadening  is  apparently  caused 
by  the  approach  of  daughter  singularities  created  from  an  initial  zero  in  |C|  >  1  (see 
Siegel,  Tan  veer  &  Dai  [40]). 

The  above  examples  show  that,  in  some  cases,  surface  tension  causes  a  singular 
perturbation  in  0(1)  time,  whereas  in  other  cases  it  does  not.  Our  method  gives 
us  a  means  of  discerning  these  cases  through  computation  of  daughter  singularity 
trajectories  in  the  extended  domain  [43].  We  remark  that  our  attempts  to  calculate 
the  5  =  0  solutions  using  the  boundary  integral  method  with  filtering  failed  before 
the  interface  had  advanced  very  far.  This  is  because  the  close  presence  of  strong 
singularities  causes  fast  growth  of  the  high  wavenumber  modes,  and  numerical  noise 
quickly  contaminates  the  computation.  Thus,  boundary  integral  methods  often  appear 
unsuitable  for  comparing  5  =  0  solutions  to  those  for  5  >  0  over  times  in  which  the 
interface  becomes  significantly  deformed.  Dai  &  Shelley  [8]  report  related  problems  in 
5  =  0  calculations,  as  discussed  in  the  introduction. 

We  turn  now  to  a  consideration  of  the  influence  of  additional  —4/3  singularities 
with  weak  amplitudes  during  finger  formation.  Figure  7a  illustrates  the  interface 
evolution  resulting  from  initial  data  of  the  form  (23)  with  Ji  =  2,  J  =  3  and  with 

^This  is  accomplished  by  prescribing  initial  data  in  Z(  of  the  form  2,;(C,0)  =  —2(1  —  C^/Cs  )^/(’''C)- 
Here  we  use  a  =  —3.2  and  =  (1.6, 0.0). 
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singularity  strengths  =  -4/3  for  j  ~  The  initial  singularity  positions 

were  chosen  as  Ci(^)  =  -(2(0  ==  (1.65,0,0),  and  ^3(0  =  (**“0.6, 1.83),  and  the  initial 
amplitudes  as  =  £'^((,0)  =  (-0.9, 0.0),  and  £^((,0)  =  (—0.01,0.03).  The 

calculation  uses  A^=512  points  and  a  time  step  At  ==  0.0005.  Two  of  the  singularities 
(Cl  and  C2)  reside  on  the  real  line,  and  these  have  large  enough  amplitude  to  produce 
the  wide  bulges  of  fluid  centered  at  y  =  1  and  y  =  -I  which  define  the  main  finger. 
As  the  third  singularity  approaches  the  unit  disk,  it  generates  a  small  bulge  of  fluid 
that  gives  the  appearance  of  the  formation  of  a  dimple  on  the  evolving  finger. 

The  dimple  is  clearly  stationary  in  the  laboratory  frame,  despite  the  overall  growth 
of  the  main  finger.  Such  behavior  has  been  well  documented  in  laboratory  experiments 
(  e.g.  Park  &  Homsy  [28],  Tabeling,  Zocchi,  &  Libchaber  [41])  and  in  numerical 
calculations  (e.g.  DeGregoria  and  Schwartz  [9,  10])  with  B  ^  0.  Our  results  show 
that  it  is  the  specific  nature  of  the  trajectory  of  the  singularity  C3(0  accounts 
for  this  behavior.  We  show  the  trajectories  of  all  the  singularities  in  Figure  7b.  At 
first,  a  dimple  starts  to  form  as  (3(0  approaches  the  unit  disk.  As  (3(0  begins  to 
move  around  the  boundary  of  the  unit  disk  towards  C  =  “1,  the  dimple  continues  to 
grows  but  remains  stationary  in  the  laboratory  frame.  With  the  assumption  that  the 
singularity  is  close  to  either  C  =  ±1,  we  can  show  that  its  speed  is  actually  the  correct 
speed  for  it  to  remain  stationary  in  the  laboratory  frame.  Details  will  be  provided 
elsewhere.  The  point  to  be  made  here  is  only  that  certain  physical  properties  can 
be  understood  in  terms  of  the  motion  of  singularities  in  the  complex  (  plane.  Under 
the  presence  of  surface  tension,  the  narrow  finger  will  eventually  widen;  however,  the 
formation  of  the  dimple  and  its  relationship  to  the  motion  of  the  singularities  will  not 
be  affected. 

Placing  additional  singularities  in  the  complex  plane  produces  additional  dimples 
and  leads  to  the  appearance  of  side-branching.  An  example  of  side-branching  due  to 
multiple  pole  singularities  is  given  in  Figure  8a,  We  use  pole  singularities,  rather  than 
a  =  -4/3  singularities,  since  our  ability  to  track  them  arbitrarily  close  to  the  unit 
disc  leads  to  a  more  dramatic  example  of  sidebranching,  but  in  reality  side-branching 
is  more  likely  to  occur  from  the  patterns  of  -4/3  singularities  created  by  the  initial 
transformation  of  zeros  in  Initial  data  corresponding  to  (27)  with  —  2  and  J  =  14 
is  used  to  generate  the  profiles.  The  starting  amplitudes  £’*^((',0)  and  singularity 
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positions  (j(0)  are  given  in  Table  I.  The  two  main  singularities  Ci  and  (2  produce  the 
finger  centered  in  the  channel,  and  remaining  singularities  (j  for  j  -  3, J  of  smaller 
amplitude  are  located  to  cause  side-branches  to  form  near  the  base  of  the  finger. 
Additional  singularities  can  be  placed  to  allow  side-branching  to  run  the  length  of  the 
finger.  As  expected,  the  dimples  and  the  corresponding  side-branches  are  stationary 
in  the  laboratory  frame.  The  motion  of  the  singularities  is  shown  in  Figure  8b,  and 
clearly  the  singularities  with  smaller  amplitudes  are  attracted  to  the  points  C  =  ±1. 
When  they  get  close  enough  to  either  of  these  points,  an  unperturbed  finger  continues 
to  grow  while  the  dimples  remain  near  its  base.  If  additional  singularities  are  present 
that  start  much  further  away  from  the  unit  disk,  they  will  arrive  close  to  the  unit  disk 
at  later  times,  producing  new  dimples  near  the  tip  of  the  finger.  The  endless  presence 
of  singularities  streaming  in  from  infinity  can  generate  fingers  with  side-branching 
patterns  seen  frequently  in  experiments. 

In  Figure  7a,  the  small  amplitude  of  the  singularity  at  (3  causes  the  formation 
of  a  small  dimple  on  the  side  of  a  well-developed  finger.  In  contrast,  when  Re 
is  the  same  order  as  Re  and  Re  the  close  approach  of  (3  to  |  (  |=  1  leads 
to  a  ‘tip-splitting’  event  and  the  formation  of  two  fingers  which  eventually  compete. 
An  example  is  given  in  Figure  9a  for  initial  data  of  the  form  (23)  with  Ji  =  2,  J  = 
3  and  Oj  =  -4/3  for  j  =  1,...,3.  The  initial  singularity  positions  are  chosen  as 
(1(0)  =  -(2(0)  =  (2.5, 0.0),  (3(0)  =  (0.0, 4.0)  and  the  initial  amplitudes  as  £^i(C,0)  = 
(-1.0, 0.0),  £2(C,0)  =  (-0.8, 0.0),  F^3(C,0)  =  (-0.3, -0.16).  We  use  N  =  512  and 
At  =  0.0005.  As  shown  in  Figure  9b,  the  singularities  initially  move  radially  towards 
the  unit  circle.  As  the  singularity  at  (3  approaches  |C|  =  1,  a  dimple  forms  near  the 
finger  tip.  As  time  advances,  Cs  gets  closer  to  the  interface  and  the  dimple  elongates 
into  a  large  indentation,  giving  the  appearance  of  two  fingers.  Eventually,  the  motion 
of  the  singularity  is  predominantly  tangent  to  the  circle,  and  the  tangential  velocity 
of  Cs  is  such  that  the  indentation  is  fixed  in  the  laboratory  frame.  The  direction  the 
singularity  moves  around  the  circle  determines  which  finger  will  dominate.  In  this  case 
the  singularity  at  (3  moves  towards  (  =  -1.  Since  the  point  C  =  “1  corresponds  to 
the  bottom  end  of  the  interface,  this  motion  has  the  effect  of  stretching  the  interface 
so  that  the  indentation  separating  the  fingers  lies  closer  to  the  bottom  end.  The  only 
way  this  can  be  done  while  also  keeping  the  indentation  fixed  in  the  laboratory  frame 
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is  for  the  upper  finger  to  grow  significantly  more  than  the  lower  one,  as  indicated  in 
Figure  9a.  Unfortunately,  the  simulation  cannot  be  not  run  long  enough  to  produce 
a  clear  outcome  in  the  finger  competition,  because  of  the  presence  of  the  singularity 
in  Ei,  and  G-  at  Nevertheless,  we  are  able  to  compute  long  enough  to  make 

the  trend  in  singularity  motion  clear. 

We  conclude  this  section  with  a  brief  examination  of  some  scenarios  which  may  lead 
to  a  self-intersection  of  the  interface.  A  self-intersection  event  is  often  referred  to  as  a 
topological  singularity;  the  possible  formation  of  this  type  of  singularity  in  Hele-Shaw 
flow  and  in  other  free  surface  flows  is  a  topic  of  much  current  interest.  A  topological 
singularity  occurs  when  the  conformal  map  ceases  to  be  univalent,  i.  e.,  when 

two  points  on  the  (  semi  circle  map  to  a  single  point.  When  this  happens,  either  the 
more  viscous  or  less  viscous  fluid  region  is  divided  into  two  disjoint  sections.  Bertozzi, 
Brenner,  Dupont,  &  Kadanoff  [7]  and  Goldstein,  Pesci,  &  Shelley  [12]  have  investigated 
possible  topological  singularity  formation  in  Hele-Shaw  flow  with  surface  tension.  They 
consider  a  particular  geometry,  consisting  of  a  vertical  ceU  with  a  thin  layer  of  fluid 
resting  on  the  bottom,  chosen  so  that  a  variant  of  the  lubrication  approximation  can 
be  applied.  Using  this  approximation,  they  concluded  that  in  certain  cases  the  top 
interface  of  the  layer  touches  the  bottom  of  the  ceU  in  finite  time.  However,  their 
geometry  and  pressure  conditions  are  significantly  different  from  ours  and  it  is  unclear 
if  their  results  can  be  extrapolated  to  our  geometry,  where  there  is  a  constant  pressure 
gradient  far  ahead  of  the  finger. 

In  some  situations,  a  loss  of  univalence  is  possible  even  when  the  singularities  and 
zeroes  of  I  C  l>  1  remain  a  finite  distance  from  |  (  |—  1.  We  used  our  code  to 

search  for  such  an  event  in  the  zero  surface  tension  problem.  Unlike  commonly  used 
boundary  integral  methods  which  run  into  resolution  difficulties  when  the  interface 
is  about  to  pinch  (see  Baker  &  Shelley  [3]),  our  numerical  approach  based  on  the 
conformal  mapping  function  will  not  incur  difficulties  as  long  as  the  pinching  is  not 
accompanied  by  a  singularity  or  zero  of  Z(  impinging  on  |  C  |=  1- 

Unfortunately,  we  were  unable  to  find  any  occurance  of  this  type  of  singularity 
in  the  situation  where  a  fluid  of  zero  viscosity  displaces  a  viscous  fluid.  Of  course, 
the  formation  of  a  topological  singularity  in  this  case  cannot  be  completely  ruled  out 
from  our  limited  examination  and  further  study  is  required.  When  we  reversed  the 
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pressure  gradient  at  infinity,  however,  so  that  the  more  viscous  fluid  on  the  right  of 
Figure  1  displaces  the  fluid  of  negligible  viscosity  on  the  left,  topological  singularities 
were  sometimes  observed.  It  is  well  known  from  the  Saffman  &  Taylor  [36]  analysis 
that  a  planar  interface  in  this  situation  is  stable  and  that  any  small  deformation  will 
reduce  with  time.  This  can  be  expected  to  be  true  even  for  most  finite  amplitude 
disturbances.  However,  our  findings  show  that  if  the  interface  is  highly  deformed 
initially,  pinching  can  occur. 

One  such  example  is  presented  in  Figure  10.  This  figure  shows  interfacial  profiles 
before  and  after  topological  singularity  formation  for  an  initial  value  problem  with 
two  branch  point  singularities  initially  close  to  the  unit  disk.  We  picked  data  in 
of  the  form  2<(Ci0)  =  -2(1  -  C^/Cs)“/(^C)  where  a  =  -1.9  and  =  (1.105,0.0); 
with  this  choice,  zeros  in  are  placed  initially  at  infinity.!!  The  change  to  liquid 
pushing  air  reverses  the  direction  of  the  characteristic  velocities  (given  by  qi),and 
the  singularities  move  outward  from  the  unit  circle.  Consequently,  the  evolution  of 
the  interface  can  be  obtained  without  explicitly  representing  the  singularities.  We 
therefore  set  2  -  (2/)r)lnC  =  G+  +  Go  and  solve  (39)  on  the  circle  |C|  =  there  is 
no  need  to  evolve  the  radius  of  this  circle.  In  the  run  we  set  N  =  128  and  At  =  0.001. 

Our  example  shows  that  a  loss  of  univalency  can  occur  in  zero  surface  tension 
flows  without  concurrent  singularities  in  z^,  and  illustrates  the  ability  of  our  numer¬ 
ical  method  to  contend  with  such  self  intersections.  The  occurance  of  these  singu¬ 
larities  appears  to  be  quite  sensitive  to  initial  conditions,  with  fatter  initial  fingers 
typically  evolving  to  a  flat  sheet  without  pinching.  The  form  of  the  solution  after  a 
self-intersection  remains  an  open  question. 

6  Discussion  and  Conclusions 

We  have  described  a  niimerical  method  designed  to  track  singularities  present  in  the 
conformal  map  from  the  unit  semi-circle  to  the  physical  domain  of  Hele-Shaw  flow  in  a 
channel.  The  method  is  restricted  to  those  conformal  maps  that  contain  singularities 
of  the  form  (23),  and  their  initial  location  Cj(0)  power  aj  must  be  given  if  they  lie 

I* Thus,  the  presence  of  surface  tension  will  not  significantly  affect  the  evolution  of  the  interface  for 
the  times  shown. 
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inside  the  circle  of  radius  i?(0).  Furthermore,  their  amplitude  E^{(,  0)  and  6'(Ci  0)  must 
be  specified,  or  at  least  their  inner  and  outer  components  must  be  given  on  |C|  =  R{0) 
and  l^l  =  1  respectively.  This  detailed  information  of  the  initial  properties  of  the 
conformal  map  correspond  to  high  precision  knowledge  of  the  initial  interface  location. 
Our  numerical  method  then  advances  the  conformal  map,  and  hence  the  interface, 
numerically  in  a  stable  way.  In  other  words,  round-off  errors  do  not  contaminate  the 
high  precision  specification  of  the  interface  location. 

Besides  the  restrictions  on  the  initial  conditions,  there  are  two  other  limitations  on 
our  method.  Singularities  in  inner  analytic  components  of  and  occur 

inside  the  unit  disk  at  As  Cj{t)  approaches  the  unit  disk,  these  singularities  also 

approach  the  boundary  of  the  unit  disk,  and  cause  slow  decay  of  the  Fourier  modes  for 
E^{^,t)  and  G{Q,t).  Consequently,  we  need  many  Fourier  modes  to  ensure  reasonable 
accuracy.  Also,  at  present  the  method  is  limited  to  5  =  0.  When  i?  7^  0,  zeros  in  2^ 
are  transformed  into  patterns  of  —4/3  singularities.  Some  of  these  can  move  toward 
the  unit  circle  and  affect  the  shape  of  the  interface  at  later  times.  Nevertheless,  as  long 
as  aj  <  -4/3  and  none  of  the  singularities  formed  out  of  initial  zeros  approach  the 
unit  disk  closely,  our  results  will  be  the  correct  limiting  behavior  as  5  0.  Perhaps 

more  importantly,  our  method  enables  accurate  B  —  0  computations  to  be  obtained 
for  quite  a  general  distributions  of  initial  singularities,  so  that  comparisons  with  B  >  0 
solutions  can  be  made.  These  kind  of  comparisons  complement  the  asymptotic  theory, 
and  facilitate  an  understanding  of  the  influence  of  small  capillary  effects. 

Despite  the  limitations,  we  find  singularities  induce  interfacial  structure  that  is 
typical  of  experiment  observations  when  B  is  very  small.  In  particular,  two  singulari¬ 
ties,  placed  on  the  real  axis  on  either  side  of  the  unit  circle,  induce  formation  of  a  long 
finger.  Singularities  off  the  real  axis  induce  small  indentations  on  this  finger  if  their 
amplitudes  are  small,  giving  the  appearance  of  side-branches,  or  large  indentations  if 
their  amplitudes  are  comparable  to  the  ones  on  the  real  axis,  giving  the  appearance 
of  tip-splitting  and  finger  competition.  In  general,  we  expect  that  a  continual  inward 
stream  of  singularities  of  all  amplitudes  can  account  for  multiple  branching  and  com¬ 
petition  as  observed  experimentally.  Although  some  aspects  of  interfacial  evolution 
due  to  multiple  singularities  have  been  examined  previously  by  Howison  [19]  using  a 
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class  of  exact  solutions  for  simple  pole  singularities,  our  algorithm  allows  a  broader 
study  to  be  undertaken  for  collections  of  isolated  singularities  of  more  general  form. 

We  plan  to  continue  studies  of  the  properties  of  the  singularities  and  the  interfa¬ 
cial  structure  their  induce.  For  example,  we  wish  to  explore  what  role  the  complex 
amplitudes  have  on  the  trajectories  of  the  singularities.  We  hope  to  find  ways  to  rep¬ 
resent  the  singularities  in  better  forms  that  may  remove  some  of  the  limitations  in  our 
method,  and  we  hope  to  find  ways  to  capture  the  transformation  of  zeros  when  B  is 
very  small,  but  non-zero.  The  authors  would  like  to  acknowledge  support  from  NASA 
grant  NAG  3-1415  (G.B.),  Department  of  Energy  contract  DE-FG02-92ER14270  (M. 
S.  and  S.  T.),  and  an  NSF  Postdoctoral  Fellowship  (M.  S.).  S.  T.  was  partially  sup¬ 
ported  by  NASA  grant  NASl-18605  while  in  residence  at  The  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE). 
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Less  Viscous  Fluid 


More  Viscous  Fluid 
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y=-i 


1.  The  unit  semi-circle  in  the  ^  plane  is  mapped  into  the  viscous  fluid  region  of  the 
channel,  with  the  circular  arc  being  mapped  to  the  interface.  The  points  A,n, 
and  C  in  the  (  plane  are  mapped  to  the  corresponding  points  in  the  channel. 
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.  A  finger  produced  by  two  branch  point  singularities  of  strength  a  =  -4/5, 
located  at  symmetric  positions  on  the  real  line.  The  profiles  show  the  position 
of  the  interface  from  i  =  0  to  <  =  0.8  in  increments  of  .1,  and  a,t  t  =  0.85.  In  the 
last  profile  (i  =  0.85)  the  singualrities  are  located  at  =  -(2  =  1.0081. 
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.  A  finger  produced  by  two  branch  point  singularities  of  strength  a  =  -4/3, 
located  symmetrically  on  the  real  line.  The  profiles  show  the  position  of  the 
interface  from  t  =  0  to  t  =  1.1  in  increments  of  .1.  In  the  last  profile  (t=l.l)  the 
singularities  are  located  at  =  —(,‘2  =  1.0054. 
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from  /,  =  0  to  /  =  1.1  in  increments  of  .1.  Only  odd  values  of  k  are  plotted, 
since  the  values  of  E^{k)  for  even  k  are  zero  by  symmetry  (in  the  computed 
spectrum  these  values  remain  at  round-off  levels).  The  values  of  E^(k)  for  fc  >  1 
are  identically  zero. 


N  0 

E 


-0.5 


-1 


-1.5 


-2  - -  ■  I - 1 - 1 _ 1  ■ 

-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 

Re  zeta 


b)  The  singularity  trajectories  corresponding  to  6(a).  After  reaching  tlie  unit 
circle,  the  (upper)  complex  singularity  Ca  moves  in  a  counterclockwise  direction. 
In  the  last  profile  (/  =  .22),  the  singularities  are  located  at  Ci  =  1.1522,(2  = 
-1.1542,  and  (3  =  (-.9381,  .3546). 
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1.2  in  increments  o: 


(b)  The  singularities  trajectories  corresponding  to  7(a). 
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-3-2-10  1  2  3 

Re  zeta 


b)  The  singularity  trajectories  for  the  finger  competition  shown  in  8(a).  The 
motion  of  the  singularities  is  towards  (  =  -1. 


47 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No  0704-0188 


Public  reporting  burden  for  t^his  collection  of  information  is  estimated  to  average  1  hour  per  response,  includingthe  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204,  Arlington,  vA  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  OHVi (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

_ February  1995 _  Contractor  Report 


4.  TITLE  AND  SUBTITLE 

A  WELL-POSED  NUMERICAL  METHOD  TO  TRACK  ISO¬ 
LATED  CONFORMAL  MAP  SINGULARITIES  IN  HELE-SHAW 
FLOW 


6.  AUTHOR(S) 

Gregory  Baker 
Michael  Siegel 
Saleh  Tan  veer 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


5.  FUNDING  NUMBERS 


C  NASI-19480 
WU  505-90-52-01 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  95-7 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

NASA  CR-195035 
ICASE  Report  No.  95-7 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor;  Dennis  M.  Bushnell 
Final  Report 

Submitted  to  Journal  of  Computational  Physics 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

U^  ncIassified-Unlimited 


12b.  DISTRIBUTION  CODE 


Subject  Category  64 


13.  ABSTRACT  (Maximum  200  words) 

We  present  a  new  numerical  method  for  calculating  an  evolving  2-D  Hele-Shaw  interface  when  surface  tension  effects 
are  neglected.  In  the  case  where  the  flow  is  directed  from  the  less  viscous  fluid  into  the  more  viscous  fluid,  the  motion 
of  the  interface  is  ill-posedj  small  deviations  in  the  initial  condition  will  produce  significant  changes  in  the  ensuing 
motion.  This  situation  is  disastrous  for  numerical  computation,  as  small  round-off  errors  can  quickly  lead  to  large 
inaccuracies  in  the  computed  solution.  Our  method  of  computation  is  most  easily  formulated  using  a  conformal  map 
from  the  fluid  domain  into  a  unit  disk.  The  method  relies  on  analytically  continuing  the  initial  data  and  equations 
of  motion  into  the  region  exterior  to  the  disk,  where  the  evolution  problem  becomes  well-posed.  The  equations  are 
then  numerically  solved  in  the  extended  domain.  The  presence  of  singularities  in  the  conformal  map  outside  of  the 
disk  introduces  specific  structures  along  the  fluid  interface.  Our  method  can  explicitly  track  the  location  of  isolated 
pole  and  branch  point  singularities,  allowing  us  to  draw  connections  between  the  development  of  interfacial  patterns 
and  the  motion  of  singularities  as  they  approach  the  unit  disk.  In  particular,  we  are  able  to  relate  physical  features 
such  as  finger  shape,  side-branch  formation,  and  competition  between  fingers  to  the  nature  and  location  of  the 
singularities.  The  usefulness  of  this  method  in  studying  the  formation  of  topological  singularities  (self-intersections 
of  the  interface)  is  also  pointed  out. 


14.  SUBJECT  TERMS 

Hele-Shaw  instability;  ill- posed;  singularities;  side-branching;  finger  competition 


15.  NUMBER  OF  PAGES 

50 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


NSN  7540-01-280-5500 


16.  PRICE  CODE 

A03 


18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 
OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

Unclassified 


Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 
298-102 


